Get some US Census data

Download US Census data on natural resource dependent jobs (total # and median income) within counties and summarize across NEP boundaries:

nep_sf <- st_read(dsn = "./data-raw", layer = "NEP_Boundaries10162018")
## Reading layer `NEP_Boundaries10162018' from data source `D:\rprojects\anep_bea_stats\data-raw' using driver `ESRI Shapefile'
## Simple feature collection with 28 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -124.7411 ymin: 18.31383 xmax: -65.87813 ymax: 48.99944
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
            "ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
            "MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
            "NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
            "SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY",
            "PR")

metrics <- load_variables(2017, "acs5", cache = TRUE) %>% 
       filter(grepl("Fishing|fishing", label))

totaljobs_sf <- get_acs(geography = "county", variables = c(Ag_For_Fish_Hunt_Min = "C24050_002"), 
                         state = states, geometry = TRUE)

nr_medinc_sf <- get_acs(geography = "county", variables = c(Med_Income = "B24031_003"), 
                         state = states, geometry = TRUE)

ustotaljobs_sf <- st_transform(totaljobs_sf, st_crs(nep_sf))
usnrmedinc_sf <- st_transform(nr_medinc_sf, st_crs(nep_sf))

nep_jobs_intersects <- st_intersects(nep_sf, ustotaljobs_sf)
nep_medinc_intersects <- st_intersects(nep_sf, nr_medinc_sf)

nep_sel_sf <- ustotaljobs_sf[unlist(nep_jobs_intersects),]
nep_sel2_sf <- nr_medinc_sf[unlist(nep_medinc_intersects),]

nep_jobs <- st_join(nep_sf, nep_sel_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)
nep_nrmedinc <- st_join(nep_sf, nep_sel2_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)

nep_jobs_sum <- nep_jobs %>% 
               select(NEP_NAME, estimate) %>%
               group_by(NEP_NAME) %>% 
               summarise(jobs = sum(estimate))

nep_medinc_sum <- nep_nrmedinc %>% 
               select(NEP_NAME, estimate) %>%
               group_by(NEP_NAME) %>% 
               summarise(medinc = median(estimate))

Plot US Census data by NEP Project Areas

Plot US Census natural resource dependent jobs data within the NEP boundaries:

pal <- colorNumeric(palette = "viridis", domain = nep_jobs_sum$jobs, n = 10)

nep_jobs_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Jobs = ",jobs),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal(jobs)) %>%
    addLegend("bottomright", 
              pal = pal, 
              values = ~ jobs,
              title = "Natural Resource Dependent Jobs</br>(Ag.,Forest.,Fishing,Hunting,Mining)",
              opacity = 1)

Plot US Census natural resource dependent jobs, median income data within the NEP boundaries:

pal2 <- colorNumeric(palette = "viridis", domain = nep_medinc_sum$medinc, n = 10)

nep_medinc_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Median Income = ",medinc),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal2(medinc)) %>%
    addLegend("bottomright", 
              pal = pal2, 
              values = ~ medinc,
              title = "Natural Resource Dependent Jobs</br>(Median Annual Income)",
              opacity = 1)

Get some US BEA data

Import US BEA data on natural resource dependent jobs (total county income) and summarize across NEP boundaries:

nr_income <- read.csv("./data-raw/CAINC5N__ALL_AREAS_2001_2017.csv", header = TRUE) %>% 
             filter(LineCode==100,str_detect(GeoFIPS, "....0")==FALSE) %>%    
             mutate(GEOID = as.character(GeoFIPS),
                    YR2017 = as.character(X2017)) %>% 
             select(GeoName, GEOID, YR2017)
nr_income$GEOID <- gsub(" ", "", nr_income$GEOID, fixed = TRUE)
nr_income$Y2017 <- as.numeric(nr_income$YR2017)

nr_inc_sf <- left_join(nep_sel_sf, nr_income, by = c('GEOID' = 'GEOID'))

nep_totinc <- st_join(nep_sf, nr_inc_sf, join = st_intersects) %>%
                   sf::st_buffer(dist = 0)

nep_inc_sum <- nep_totinc %>% 
               select(NEP_NAME, Y2017) %>%
               group_by(NEP_NAME) %>% 
               summarise(totinc = sum(Y2017, na.rm = TRUE))
nep_inc_sum$totinc <- ifelse(nep_inc_sum$totinc==0,NA,nep_inc_sum$totinc*1000)

Plot US BEA natural resource dependent jobs total income estimates within the NEP boundaries (where available):

pal3 <- colorNumeric(palette = "viridis", domain = nep_inc_sum$totinc, n = 10)

nep_inc_sum %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Positron") %>%
    addPolygons(popup = ~ paste(NEP_NAME,"</br>","Total Income (2017) = ",totinc),
                stroke = FALSE,
                smoothFactor = 0,
                fillOpacity = 0.7,
                color = ~ pal3(totinc)) %>%
    addLegend("bottomright", 
              pal = pal3, 
              values = ~ totinc,
              title = "Natural Resource Dependent Jobs</br>(Total Annual Income, if reported)",
              opacity = 1)